Maximally Localized Wannier Functions within the FLAPW formalism 
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We report on the implementation of the Wannier Functions (WFs) formalism within the full- 
potential linearized augmented plane wave method (FLAPW), suitable for bulk, film and one- 
dimensional geometries. The details of the implementation, as well as results for the metallic SrVC>3, 
ferroelectric BaTi03 grown on SrTiC>3, covalently bonded graphene and a one-dimensional Pt-chain 
are given. We discuss the effect of spin-orbit coupling on the Wannier Functions for the cases of 
SrVC>3 and platinum. The dependency of the WFs on the choice of the localized trial orbitals as 
well as the difference between the maximally localized and "first-guess" WFs are discussed. Our 
results on SrVC>3 and BaTiOs, e.g. the ferroelectric polarization of BaTiOs, are compared to results 
published elsewhere and found to be in excellent agreement. 
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I. INTRODUCTION 

Commonly, the electronic structure of periodic solids 
is described in terms of Bloch functions (BFs), which are 
eigenf unctions of both the Hamiltonian and lattice trans- 
lation operators. Due to their delocalized nature BFs are 
difficult to visualize and hence do not offer a very intu- 
itive picture of the underlying physics. Furthermore, BFs 
do not provide an efficient framework for the study of lo- 
cal correlations. An alternative approach to electronic 
structure that does not exhibit these weaknesses is pro- 
vided by maximally localized Wannier functions (ML- 
WFs). Related to the BFs via a unitary transformation, 
MLWFs constitute a mathematically equivalent concept 
for the study of electronic structure. They are well lo- 
calized in real space and in contrast to the complex BFs 
purely real^. Therefore, it is easy to visualize them and 
to gain physical insight e.g. into the bonding properties 
of the system under study by extracting characteristic 
parameters such as the MLWFs' centers, spreads, and 
hopping integrals as well as by analyzing their shapes. 

Wannier functions (WFs) were first introduced by 
Wannier in 1937- as the Fourier transforms of BFs. Sim- 
ilar to a ^-function, which is the Fourier transform of 
a plane wave, WFs are localized in real space while the 
BFs are not. However, BFs are only determined up to an 
arbitrary phase factor, and hence the definition of WFs 
as Fourier transforms of BFs does not specify the WFs 
uniquely. As the localization properties of the WFs de- 
pend strongly on the phase factors of the BFs, the Wan- 
nier function approach experienced little enthusiasm un- 
til very recently, after methods for the calculation of WFs 
with optimal localization properties had been developed. 
One of these new techniques for the construction of lo- 
calized WFs is based on the N-th order muffin-tin-orbital 
(NMTO) method^! 4 -! 5 - Another method performs at each 
fe-point a unitary transformation among the BFs belong- 
ing to different bands yielding a new set of functions, 
the Fourier transforms of which are the MLWFs. 6 The 
MLWFs approach is not limited to insulators but is also 



capable of providing well localized orbitals for metalsi 7 - 
Only the latter technique is considered in this work. 

Sheding new light on otherwise hard to calculate 
properties of materials, nowadays MLWFs have almost 
reached the popularity of BFs, and using both allows 
to achieve a rich diversity in understanding, originat- 
ing from revealing both itinerant and localized aspects of 
electrons in periodic potentials. For example, a modern 
theory of polarization 8 -" ^ 10 ' 11 ' 12 is based on the displace- 
ments of the centers of the MLWFs. The orbital polariza- 
tion may be expressed in terms of MLWFs i 13 ' 14 Studying 
the MLWFs for disordered systems yields a transparent 
description of bonding properties. 15 MLWFs provide a 
minimal basis set that allows for efficient computations 
of the quantum transport of electrons through nanos- 
tructures and molecules. 16,17 Within the research area of 
strongly correlated electrons MLWFs are becoming the 
preferred basis for studying the local correlations , 1 8 1 1 9 1 20 

The MLWFs-induced burst in studying the properties 
of materials which are hard to probe on the basis of 
traditional band theory is very recent and many sub- 
tle aspects, such as magnetism, various spin-orbit cou- 
pling and non-collinearity-driven effects are still to be 
put on the MLWFs footing. In this respect the precision 
of the computational electronic structure method used 
for the construction of the MLWFs might play a very 
important role, as sophisticated details of the electronic 
structure and tiny energy scales are involved. In particu- 
lar in magnetism, the choice of the appropriate ab initio 
method plays a crucial role. From this point of view it is 
common consensus that the full-potential linearized aug- 
mented plane wave method (FLAPW) is one of the most 
precise electronic structure methods used today. Ab- 
initio MLWFs have already been calculated within the 
FLAPW framework for MnO 21 and Ti0 2 M£l 

In the present paper we report in detail on the im- 
plementation of MLWFs within the FLAPW method as 
implemented in the FLEUB 2 ^ code. The current imple- 
mentation allows a fast computation of MLWFs for a 
large variety of materials and complex geometries, in- 
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eluding bulk, film 25 and truly onc-dimcnsional geomet- 
rical setups^ To verify our implementation we apply 
the method to four different systems, two different per- 
ovskite systems, SrVC>3 and BaTiC>3, one metallic and 
one ferroelectric, graphene, a covalently bonded material, 
and a one-dimensional Pt-chain. This article is struc- 
tured as follows: We start in section II with a short out- 
line of MLWFs and their construction procedure, defin- 
ing the quantities required from the first-principles cal- 
culation based on the density functional theory (DFT) 
by the maximal localization algorithm. First-guess WFs 
- originally devised as a starting point for the MLWF- 
algorithm, but widely used as a suitable alternative to 
the MLWFs - are introduced. Then, the details of our 
FLAPW implementation are described. In Section III 
we apply the formalism to SrVC>3, BaTi03, graphene 
and a one-dimensional Pt-chain. We discuss the effects 
of spin-orbit coupling on the MLWFs for SrVC>3 and the 
Pt-chain. We compare our results on SrVC>3 and BaTi03 
with theoretical and experimental data, respectively, and 
find excellent agreement. Finally we close this work with 
conclusions in Section IV. 



II. METHOD 

A. Maximally localized Wannier functions 

For an isolated band, i.e. a band that does not become 
degenerate with other bands at any fc-point, with cor- 
responding BFs \tffk), the definition of WFs as Fourier 
transforms of BFs leads to the following expression: 

k 

where R is a direct lattice vector, which specifies the 
unit cell the WF belongs to, and the Brillouin zone is 
represented by a uniform mesh of N fc-points. The |^k) 
are normalized with respect to the unit cell, while the 
\Wr) constitute an orthonormal basis set with respect 
to the volume of N unit cells. 

However, Eq. (JTj) does not define the WFs uniquely: 
The BFs are determined only up to a phase factor - 
hence, for a given set of BFs and a general fc-point de- 
pendent phase 0(k), 

I^R>' = ^E e ~ ik ' Rei<Kk) l^> ( 2 ) 

k 

equally constitute a set of WFs. For their use in prac- 
tice, it is desirable to have WFs that decay exponen- 
tially in real space, exhibit the symmetry properties of 
the system studied, and are real- rather than complex- 
valued functions^. For the one-dimensional Schrodinger 
equation and an isolated single energy band, Kohn2£ has 
shown that there exists only one WF which is reali, falls 
off exponentially with distance and has maximal symme- 
try. WFs with maximal spatial localization 6 (MLWFs) 



fulfill these requirements of real-valuedness 1 , optimal de- 
cay properties and maximal symmetry. The constraint 
of maximal localization eliminates the nonuniqueness of 
WFs and determines 0(k) up to a constant. 

In the general case, energy bands cross or are degener- 
ate at certain fc-points, making it necessary to consider 
a group of bands. This increases the freedom in defining 
WFs further, as now bands may be mixed at each fc-point 
via the transformation Umn'. 

\W Rn ) = l£ e - lk - R 5>i k )|^ km }, (3) 

k m 

where the BF has a band index m, the WF an orbital in- 
dex n, and the number of bands - which may depend on 
the fc-point - has to be larger than or equal to the num- 
ber of WFs that are supposed to be extracted. Impos- 
ing the constraint of maximal spatial localization on the 
WFs determines the set of {/^-matrices up to a common 
global phase In case the number of bands is equal to 
the number of WFs, the Umn matrices are unitary. This 
situation usually occurs when an isolated group of bands 
may efficiently be chosen for the system under study. In 
the more general case of entangled energy bands, 7 how- 
ever, the number of bands is fc-point dependent and Umn 
no longer unitary. 



B. Maximal localization procedure 

Requiring the spread of the WFs to be minimal im- 
poses the constraint of maximal spatial localization. For 
the spread of the WFs the sum of the second moments, 

n = £[<x 2 ) n -«x) n ) 2 ], (4) 

a 

is used, where (}„ denotes the expectation value with re- 
spect to the Wannier orbital |Won) an d the sum includes 
all WFs formed from the composite group of bands. Min- 
imization of the spread yields the set of optimal Umn- 
matrices. 

An efficient algorithm for the minimization of the 
spread Eq. has been given by Marzari and Vanderbilt 
first for isolated groups of bands, 6 and later on general- 
ized for the case of entangled energy bands The corre- 
sponding computer code is publicly available 2 ^ and was 
used in this work. Two quantities are required as input 
by this computational method and have to be provided 
by the first-principles calculation: First, the projections 
Ami. = (ipkm\gn) of localized orbitals \g n ) onto the BFs 
are needed to construct a starting point for the iterative 
optimization of the MLWFs. Second, the overlaps be- 
tween the lattice periodic parts Uk m (x) — e~ lk ' x i/>km(x) 
of the BFs at nearest-neighbor fc-points k and k + b, 

= (ukm|«k+b,n), are necessary to evaluate the 



3 



relevant observables 6 : 

< x )« = -^E^b31nM^ b ) (5) 

k,b 

and 

< x2 )« = Jf E ™ b t 1 " l^- b) | 2 + (9kM& b >) a ], (6) 

k,b 

where % is a weight associated with b, and 

M( k . b ) = )*[/(k+b) M (k,b) f7 N 

mi 7712 

evolves during the minimization process due to the itera- 
tive refinement of the Umn- The relations Eqns. are 
valid for uniform fc-point grids, while in the continuum- 
limit the fc-space expressions for the matrix elements of 
the position operator are given by£ 

(Wn n \x\W 0m ) = i ^S J rf 3 fce ik ' R (M kn |V k |M krn ) (8) 
and 

(W Rn |x 2 |W 0m ) J 'd 3 fce lk R (u k „|V^|u k „i). (9) 

Replacing the gradient V k by finite-difference expressions 
valid on a uniform fc-point mesh, one obtains the weights 
w h in Eqns. (j5j[H]). Through Eqns. (J5j[ll[7]) the spread f2 
in Eq. (j4|) may be expressed in terms of and be minimized 

(k) 

with respect to the l^Ti-matrices. 

C. First-guess Wannier functions 

The iterative optimization process requires as a start- 
ing point first guesses for the MLWFs. In order to con- 
struct these, one projects localized orbitals \g n ) onto the 
BF-subspace: 

|<M = ItfkmXVfcmbn) = E A ™« l^ k ™>- ( 10 ) 
771 rn 

As the first-guess WFs are supposed to constitute an or- 
thonormal basis set, the |0 kn ) are orthonormalized via 
the overlap matrix Smn = (<£ k m|<^ k n) 

|^k„)-^((5( k ))^)„ m |0 km ), (11) 

m 

before the WFs are calculated from them 

i^R„) = i^ e -*- R iv; kn ). (12) 

k 

While the first-guess WFs are dependent on the choice 
of localized orbitals \g n ) they converge to the one and 



only one set of MLWFs in the course of the minimization 
procedure. 

Although the first-guess WFs of Eq. (TTSj) are not 
unique they agree well with the MLWFs in many cases. 
Examples where there is substantial difference between 
first-guess WFs and MLWFs include systems where the 
centers of the Wannier orbitals do not coincide with the 
centers of the atoms. If for the system under study the 
first-guess WFs are already satisfactory, one may skip 
the localization procedure and take Eq. (fT2"]) as the final 
result. Computing WFs in such a way requires much less 
time, as the Mmn^ matrix elements do not have to be 
calculated and the minimization of the spread functional 
is not performed. First-guess WFs have been successfully 
applied to SrV0 3? ±2 Y 2 0^. and NiO^ for example. 

D. Calculation of A/^ b) within the FLAPW 
formalism 

For the calculation of MLWFs the most important 
quantity is the M^n matrix, which - according to 
Eqns. @ || - contains all information needed to de- 
termine spreads and centers. With the lattice peri- 
odic part Ukm(x) being related to its BF by u km (x) = 
e~ jk ' x ?/; km (x), tne matrix elements assume the 

form 

M<& h) - / e- lb - x (^ k „i(x))*V[ k+ b],„(x)d 3 a; . (13) 

By [k] we denote the wave vector obtained from k by 
subtracting the reciprocal lattice vector that moves k into 
the first Brillouin zone, according to [k] = k — G(k). 

Within FLAPW ) 30 i 31 space is partitioned into the 
muffin-tin (MT) spheres centered around atoms \x and 
the interstitial (INT) region. Consequently, has 
contributions from both, 

Mi k „ b) = Mi k - b ) | INT + £ Mi k ' b )| MT , , (14) 

which will be discussed separately in the following. The 
treatment of the vacuum regions occurring in film and 
one-dimensional setups is discussed in the appendices lAl 
andlBl respectively. 

Inside the muffin-tin, the BF is expanded into spherical 
harmonics, radial basis functions ui, which are solutions 
of the scalar relativistic equation at band-averaged ener- 
gies, and the energy derivatives iii of the uf. 

^km(x)|MT" 

= E(^U( k KM + Bt m (k)uf(r))Y L (r), ^ 

L 

where atom [i is located at r M and r = x — t m . Here, 
to is the band-index and L = (l,l z ) stands for the an- 
gular momentum quantum numbers I and l z . The case 
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where the lapw basis is supplemented with local orbitals 
is treated in the appendix Id Using the Rayleigh plane 
wave expansion 

e-* b ' x - 47re- iW * Y,(-l) l i l ji(rb)Y L (b)YZ{r ), (16) 



the contribution Mmn\wri' of the muffin-tin region of 
atom /i to the Mmn^ matrix reads: 

xE(( A i m ( k ))*^,n([ k + b ])*u( b ' i ' i ') 

+ (Al m (k))*^,„([k + b]K 2 (b,L,L') ll7! 
+ (B^ m (k))*^,„([k + b])^(b,L,L') 



L.L' 



+(B£ im (k))* J B£, > „([k + b])^ 2 (b,L,LO). 

The matrix elements ^(b, L",L) and t^ 2 (h, L" , L) are 
given by the sums over radial integrals 

= ^GS w '(b) | r 2 iKr&)<(r)<»(r)dr, 

*f 2 (b,i",i) (18) 
= E G W' V 'W / r 2 iK^)<(r-X(r)dr, 

L' ^ 

and analogously for t 21 and t 22 , where 

G™"(b) = G%tf m "i l '(-l) l 'Y L ,(b), (19) 
with the Gaunt coefficients 

GtfPy = / )Y l r m , {v)Y*, m „ (f ) dfi. (20) 

The quantities defined in Eq. (fl~8|) depend on the vec- 
tors b joining a given A-point to its nearest neighbors. As 
a uniform fc-mesh is used the set of b vectors and hence 
also the integrals defined in Eq. (fT8|) are independent of 
the fc-point. Thus, the quantities Eq. (fl~8|) have to be 
calculated only once. 

Employing the expansion of the BF in the interstitial 
region 



i(k+G)-x 



(21) 



V'km(x) = -j=y^c km (G)e 

V V Q 

the INT contribution to the Mmn^ matrix is deduced: 
Mi k „ b) |iNT (c k , m (G))* C[k+b] ,„(G') 



V 



G,G' 



e »([k+b]+G')-x e -i(k+G)-x e -ib x ^3^ 



(22) 



INT 



where the integration stretches over the interstitial only. 
Introducing the step function 6(x), that cuts out the 



muffin tins, and its Fourier transform 9g, Eq. ((22)) can 
be cast into the final form 



M ran I INT 



£ (c k , m (G))* C[k+b] ,„(G')e G(k+b) 

+G-G' 



(23) 



G,G' 



where G(k + b) denotes the reciprocal space vector that 
moves (k + b) into the first Brillouin zone, [k + b] = 
k + b G(k + b). 



E. 



Calculation of A%1 within the FLAPW 



formalism 



For the localized orbitals \g n ) required to determine 
the first-guess WFs, we mostly use functions that are 
zero everywhere in space except in the muffin-tin sphere 
of that atom, to which the resulting WF is attributed 
in this sense. In practice, this works not only for WFs 
that are atom-centered but also for bond-centered ones. 
Thus, g ra (x) is given by 



ff™( x ) = ^2 c nXUi{r)Y L {r), 



(24) 



where r = x — -ry, is the position relative to the center 
of the atom, to which the first-guess WF is attributed, 
and the coefficients c n ^ control the angular distribution 
of <7„(x). For the radial part ui(r) of the localized orbital 
we use the solution uf(r) of the radial scalar relativistic 
equation for the actual potential obtained from the first- 
principles calculation at an energy corresponding to the 
bands from which the WF is constructed. It is also pos- 
sible to use Gaussianspi or the radial parts of hydrogenic 
wave functions for ui{r). Where angular momentum is 
concerned in Eq. (|24p , contributions of different angular 
momenta have to be summed in the general case to allow 
the definition of hybrids such as sp 3 orbitals, while there 
is only an I = 2 contribution for WFs corresponding to d 
orbitals, for example. 

For a general radial part ui(r) the projection of the 
localized orbital \g n ) onto the BF is given by 



^ = E^[«-( k ))* / u?(r)Mr)r 2 dr 

L J 



(25) 



where the expansion of the BF given in Eq. (fTS")) was used. 
Choosing ui(r) = wf (r) Eq. (|25p simplifies as follows: 



A%1 = (VfcmlSn) = 



L 



c ™,i( a L,m( k ))* 



(26) 



In order to construct better first guesses for bond- 
centered WFs \g n ) may also be constructed as a linear 
combination of two localized orbitals - one orbital for 
each atom participating in the bond. Calculating the 
WFs for graphene in the next section we proceeded this 
way. 
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F. Wannier Representation of the Hamiltonian 

Formulating the Hamiltonian in terms of WFs is a 
particularly useful starting point when effects of correla- 
tio n 19 i 20 ' 29 are studied by DMFT. Furthermore, the hop- 
ping integrals - along with the MLWFs' spreads, centers 
and shapes - provide intuitive insight into the electronic 
structure. 

Written in terms of BFs the Hamiltonian H assumes 
the diagonal form 



be eigenstates of the projection of the spin-operator onto 
the spin-quantization axis. This means that for given n 
c n ha may differ from zero only for one spin component 
a. 

Eq. (|28p remains valid in the case of spin-orbit cou- 
pling, but the matrix elements H m ^ m i(Ri — R2) in Eq. 
([28| correspond to hopping between spinor-valued Wan- 
nier orbitals then, where the two spin-components are 
given by 



H = i e « ( k ) I ^ k ™ ) (^ k " 



I WW) = \a)(a\W Rm ), <7=U 



(32) 



(27) 



where e n (k) stand for the eigenvalues of H . If the number 
of bands is equal to the number of MLWFs extracted the 
C/rm-matrices in Eq. ([3]) are unitary. In this case we 
arrive at the equivalent form of the Hamiltonian 



Alternatively, the hopping matrix elements may be de- 
composed according to the spin-channels: 

Cm'( R l _ R 2) 

= -^^£n(k)(VF Rim(T |* k „)(* k „|WR 2m ^) 



kn 



H = E H ^nA^i - *L2)\Wn im )(Wn 2 m>\, (28) =^EE e«( k K k ' (Rl * 



(33) 



Rim R2m/ 



kn n'n' 1 



where 



x rr/ (k) ro (k) o (k) <y (k) 

\ n"m> n" no Tin' a' n' 'm' ' 



H m ,m' (R-l — R2) 



where the overlap (^knal^kn'o) is denoted 0^, a . The 
corresponding real-space representation of the Hamilto- 
nian is given by 



(29) 



N 



E^ k > 



k-(R!-R 2 ) { T r(k)\ rr(k) 



(34) 



The hopping integrals H m ^ m i(Ri — R2) quantify the 
hopping of electrons from Wannier orbital |Wr 2TO ') into 
Wannier orbital |Wr iTO ). 



G. Spin-orbit coupling 

In the case of spin-orbit coupling Eq. (flU)) assumes the 
form 

= J2 I ^ 4b - x (V-k mff (x))> [k+b] ^(x)d 3 x, (30) 

where f/'kmo-(x) is the BF with lattice vector k, band 
index n, and spin index a. The spin index a refers to the 
eigenstates of the projection of the spin-operator onto 
the spin-quantization axis. Likewise Eq. (|25p has to be 
changed into 



^ = E EE 

Rim cr,a' 

Hm m'(^-l — R-2)|W / Rim(r) (H^Ram'cr' 



Compared with Eq. (f2"9"|) the decomposition Eq. 
of the hopping matrix elements into spin-channels gives 
further insight into how the spin-channels are coupled. 

The angular characters of the spin-orbit induced cor- 
rections can be understood easily, by applying the L • S 
operator on the MLWFs that one would obtain in a cal- 
culation without spin-orbit coupling. It is convenient to 
make use of the identity 

L-S = L z S z + hL + S- + L-S + ]. (35) 



As a detailed example we consider the effect of L • S on 
T): 



L a 

xlMmAWj u» a (r)uUr)r 2 dr (31) 
+(KmAW I uL(r)u lia (r)r 2 dr]. 



L z S z \d xy }\ t> = -i\d x 2^ y 2}\ |) 



(36) 



^\d xz )\ I) - -\d yz )\ I) 



In the regime from weak to modest spin-orbit coupling 
it is reasonable to choose the localized orbitals \g n ) to 



Hence, the resulting idealized MLWF has an up- 
component the real part of which is d xy and the imag- 
inary part of which is —d x i_ y i. The real part of the 
down-component is — \d yz while the imaginary part of 
the down-component is given by —\d xz . In Table [J we 
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TABLE I: Angular part of idealized spin-orbit coupled ML- 
WFs. Columns 2,3 and 4: Components of the angular func- 
tion obtained by applying L • S to the angular function in 
column 1. 



|, real part 


T, imaginary part 


i, real part 


\, imaginary part 


dxz 

d3y2_ r 2 

Pz 


d x 2 — y 2 
2 \^3d x y 

0.0 


— 2^V Z 

d x 2 — z 2 
0.0 

\Vx 


2 dxz 
2~dxy 

— \\Fidyz 

\Vv 



list the results for various angular functions for later ref- 
erence in the results section. By 



1 



and 



-d 3z 2_ r 2 



2 d x 2 — y 2 



— V3d x 



-V3d 3z 2^ r 2 



(37) 



(38) 



we denote the angular functions obtained by rotating 
d 3z 2 _ r 2 and d x 2 _ y 2 around the x-axis by an angle of 
respectively. 

For later reference we consider the example of the Wan- 
nier orbital d xy \ T)s qa , which is an eigenstate of the pro- the three ti g bands we constructed three MLWFs, d x 



applying the generalized gradient approximation (GGA) 
to the DFT. SrV0 3 , and BaTi0 3 where calculated in 
the bulk mode of the FLEUR program, graphene in the 
film mode. For the calculation of the Pt-chain the one 
dimensional version of the program was used. 



A. SrV0 3 

The transition-metal oxide SrVC>3 crystallizes in a per- 
fectly cubic perovskite lattice with a lattice constant of 
7.26 a.u.. The Sr ions are placed at the corners of a cube 
(see Fig. [2]). The O ions are placed at the face centers 
and form an ideal octahedron in the center of which the 
V ion is located. SrVC"3 is a metal with an isolated group 
of three ti g bands around the Fermi level, which are par- 
tially occupied by one ci-electron (See Figure [1]). Within 
our GGA calculation we obtained a bandwidth of 2.5 
eV for the i2 9 -group. The experimental lattice constant 
was assumed. We used the exchange-correlation poten- 
tial of Perdew, Burke and ErnzerhofiH For Sr, V, and 
O muffin-tin radii of 2.8 a.u., 2.1 a.u. and 1.4 a.u. were 
used, respectively. Calculations were carried out with a 
plane wave cut-off of 4.5 a.u. -1 . A uniform 16x16x16 
/c-point mesh was used for the Wannier construction. For 



jection of the spin operator onto the spin-quantization 
axis. If the spin-quantization axis does not coincide with 
the z-dircction, a transformation from the states |<7) sqa 
to the basis of eigenstates of the z-component of the spin- 
operator is required before Eq. (|35p can be applied. For a 
general spin-quantization axis specified in terms of angles 
9 and <fi the transformation matrix is given by: 



d yz and d xz , 



cos I 



sin ( 



e 



sin (|) e % ~- 
— cos (!) e 1 



(39) 



After application of Eq. (|35| the states are transformed 
back to the original basis. We give the result for the 
spin-quantization axis pointing in [11 Indirection: 

LzS z d X y\ T)sqa 

— y^^ 2- ^ ^ ^ scla i , \j~~^d' X 2 —y 2 \ |)sqa 
— [L + S- + L-S+]d xy \ |)sqa 

= ^ [ d yz - d xz ] I t)sqa + [d yz +d xz } \ |) sqa 



(40) 



^~j~2~ \^ xz — d"yz\ I -L)sqa- 



III. RESULTS 

We have performed first-principles calculations within 
the framework of the density functional theory (DFT) 



which are equivalent due to symmetry. The 
MLWFs are centered at the V site. The spread, Eq. ([4J, 
of the MLWFs was found to be 6.97 a.u. 2 for each of 
the three orbitals. The first-guess WFs are characterized 
by a spread which is only 3-10 -4 a.u. 2 larger, showing 
that MLWFs and first-guess WFs are nearly identical in 
this case. To investigate the influence of spin-orbit cou- 
pling on the MLWFs a calculation including spin-orbit 
coupling was performed for the plots (see section III G[) . 
The spin-quantization axis, which defines the two spin- 
components of the spinor- valued MLWF, was chosen in 
[111] direction, to ensure that the spin components of the 
6 spin-orbit MLWFs are related by symmetry. The spin- 
orbit MLWFs are complex- valued. The imaginary parts 
of the up and down-components of the d xy \ |)-dominated 
orbital, for example, are d x 2 _ y 2 -like plus an admixture 
of d yz -d XZl while the real part of the down-component 
is (d yz + dj; 2 )-like. This result can be understood from 
the simple model in section Hi Gl that leads to Eq. (|40|) . 
The isosurface-plot for the d^-dominated orbital given in 
Fig. [2] clearly shows the hybridization between the V(i2g) 
and 0(2p) orbitals. The symmctry-inequivalent hopping 
integrals 77 m m /(Ri — R2), Eq. (j2"9")l , are listed in Ta- 
ble |TT] and found to agree well with recently published 
WF-result o 4 ' 20 on SrV03. For reasons of symmetry the 
lst-nearest-neighbor hopping integrals between different 
orbitals (e.g. d xz and d yz ) are zero in Tabic HT1 However, 
there is a coupling between the d xz orbital and the d yz 
orbitals at the 110 and 111 sites, for example. Due to 
the dominance of the nearest-neighbor hopping the three 
MLWFs may, nevertheless, approximately be considered 
independent. The fast decay of the hoppings with dis- 
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r m x r r 



FIG. 1: Bandstructure of SrV03. Red: i2 ff -bands around the 
Fermi level. 




FIG. 2: Isosurface plot of the t 2g -like MLWF d xy for SrVOs 
calculated with spin-orbit coupling. Left: Spin-up component 
(real part), isosurface=±0.05. Right: Spin-down component 
(imaginary part), isosurface=±0.001. The color of the iso- 
surface refers to the sign: Positive for dark red and negative 
for dark blue. Red balls: O sites, cyan balls: Sr sites, V 
site at the center. The WFs were plotted using the program 
XCrySDen^. 



tance furthermore indicates the short-range bonding in 
SrVC>3. The dominance of the OOT-hopping for the d xz - 
orbital over the 010-hopping reflects the restriction of 
electron hopping to the xz-plane. 

In order to study the convergence of the MLWFs with 
number of fc-points we performed a second calculation 
using an 8x8x8- mesh of /c-points. This yielded hop- 
pings identical to those of the previous calculation, but 
a slightly smaller spread of 6.73 a.u. 2 per orbital. This 
latter difference is attributed to the fact that the spread 
was calculated via the finite difference formulae Eqns. ([5j 
i). 



TABLE 


II: Hoppinj 


; Inte 


jrals 


for SrVOs. 


Ener;; 


pes are in 


meV. 














xyz 


001 


010 


011 


101 110 


111 


002 020 


dxz? dxz 


-262.0 - 


27.0 


5.8 


-84.0 5.8 


-5.7 


7.6 0.2 




0.0 


0.0 


0.0 


0.0 9.2 


3.6 


0.0 0.0 



B. BaTi0 3 

As a simple application of the Wannier-function 
scheme we present the calculation of the ferroelectric po- 
larization of the ferroelectric perovskite BaTiC>3. The 
evaluation of the polarization from a DFT calculation of 
an infinite crystal can be achieved by means of the Berry- 
phase technique. After the construction of MLWFs for 
the occupied valence bands this leads to the following 
expression for the polarization^ ' 9 ' 10 ' 11 ' 12 

P = 5>X 4 + ]T e <x} n , (41) 

i n 

where qi and Xj denote charge and position of the ion 
cores and (x)„ are the centers of the occupied Wannier 
orbitals. 

We applied this formalism to strained BaTiC>3 which is 
assumed to have been grown epitaxially on top of SrTiC>3 
assuming the in-plane lattice constant (a = 7.46 a.u.) 
of SrTiC>3. We did not consider any finite thickness or 
interface effects but simply assumed that this epitaxial 
relation will hold for reasonably thin films. The lattice 
constant perpendicular (c) as well as the positions of all 
atoms in the unit-cell where then relaxed by a series of 
force and total energy calculations. For Ba, Ti and O, 
muffin-tin radii of 2.2 a.u., 2.0 a.u. and 1.3 a.u. were 
used, respectively. The plane wave cut-off was chosen to 
be 4.8 a.u. -1 . Using the exchange correlation potential 
of Perdew and Wang2i we obtained a c/a ratio of 1.07, 
in reasonable agreement with experimental data4£ The 
resulting atomic positions are given in Table IIIII and the 
crystal structure of BaTiC-3 is illustrated in Figures[3]and 
ID Compared to the cubic perovskite structure, the oxy- 
gen atoms are moved out of the face centers and the cube 
is elongated in z-direction. z-reflection symmetry is lost. 
Az in Table IIIII specifies the displacement of the oxygen 
and titanium atoms from the symmetric positions in the 
face centers and the center of the cuboid, respectively. 

We calculated MLWFs separately for the 9 oxygen p- 
bands, the 3 barium p-bands, the 3 oxygen s-bands, the 
one barium s-band, and the 3 titanium p-bands (the re- 
maining electrons were treated as core electrons) using 
a uniform fc-point mesh of 16x16x16 fc-points. As final 
spread, Eq. |g}, 48.03 a.u. 2 were obtained for the 9 oxy- 
gen p MLWFs while the spread of the first-guess WFs 
was 48.08 a.u. 2 , demonstrating that first-guess WFs and 
MLWFs are nearly identical for BaTiC>3. Figures [3] and 0] 
show the isosurfaces of the resulting MLWFs. The ML- 
WFs clearly reflect the broken z-reflection symmetry. Ta- 
ble IIVI lists the coordinates of the centers of the MLWFs 
along with their deviations Az from the ion sites. As 
evident from there, the oxygen-MLWFs for the site close 
to the xy-plane exhibit the largest response to the bro- 
ken z-reflection symmetry. Applying Eq. (|4ip we find a 
polarization of 48.9 /xC/cm 2 in excellent agreement with 
experimental data^ of 43 /iC/cm 2 for the case of thin 
BaTiC>3 layers grown on SrTi03. The displacements of 
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TABLE III: Positions of the Ba, Ti and ions in the con- 
strained ferroelectric perovskite BaTiOs (atomic units). For 
the O ions, Az is the displacement from the face centers. For 
the Ti ion, Az specifies the displacement from the center of 
the cuboid. 





X 


y 


z 


Az 


Ba 


0.000 


0.000 


0.000 


0.000 


Ti 


3.730 


3.730 


3.901 


-0.092 


O 


3.730 


3.730 


0.449 


0.449 


O 


3.730 


0.000 


4.284 


0.292 


O 


0.000 


3.730 


4.284 


0.292 



TABLE IV: BaTiOs: Coordinates, displacements and spreads 
of the MLWFs (atomic units). 

x y z Az (x 2 ) 



4.75 
5.69 
5.69 
5.69 
5.53 
4.73 
5.69 
4.73 
5.53 
6.03 
6.15 
6.15 
2.77 
2.64 
2.64 
3.20 
1.48 
1.47 
1.47 



the centers of the MLWFs with respect to the centers of 
the atoms contribute 36% to the polarization. 

In order to assess convergence of the results with re- 
spect to the number of fc-points a comparative calcula- 
tion was performed using an 8x8x8 fc-point mesh. This 
calculation yielded a final spread of 47.19 a.u. 2 for the 
MLWFs of the 9 oxygen p bands and a total polarization 
of 48.6 /iC/cm 2 . We assume these small differences to 
be finite difference errors introduced by using formulae 
Eqns. ©El). 

C. Graphene 

Graphene is a covalently bonded system. Conse- 
quently, one expects that the MLWFs are bond centered. 
This is a particularly stringent test for our implementa- 
tion as the LAPW basis functions in which the BFs are 
expanded (see Eq. ITS)) are centered around the atoms. 
Actually, the four valence bands do not constitute an iso- 




FIG. 3: MLWFs 0(p z ) and 0(p y ) for the oxygen site close to 
xy-plane in BaTiOs. Isosurface=±0.05. Red balls in the face 
centers: O sites, cyan balls at the corners: Ba sites, green ball 
at the center: Ti site. The O site above the upper face of the 
cuboid is not depicted. 




FIG. 4: MLWFs 0(p z ), 0(p x ), and 0(p v ) for the oxygen site 
close to xz-plane in BaTiOs. Isosurface=±0.05. Red balls in 
the face centers: O sites, cyan balls at the corners: Ba sites, 
green ball at the center: Ti site. 



lated group of bands as they touch an unoccupied band 
at the JT-point. Avoiding the fsT-point when choosing 
the uniform fey-mesh, disentangling is not necessary, how- 
ever. A single layer of graphene was calculated within 
the FLEUR film mode. The muffin-tin radii and the plane 
wave cut-off were chosen to be 1.28 a.u. and 4.6 a.u. , 
respectively. The C-C bond length was assumed to be 
2.72 a.u.. We used the exchage-correlation potential of 
Perdew, Burke, and Ernzerhof.— MLWFs and first-guess 
WFs were constructed for the four valence bands using 
an 8x8 fey-mesh in the two-dimensional Brillouin zone. 
For the construction of the first-guess WFs, two calcu- 
lations were performed: In one calculation the localized 
functions \g n ) corresponding to the sp 2 -bonds were cho- 
sen to be restricted to the muffin-tin sphere of only one 
atom (FWF1), while they were restricted in the second 
calculation (FWF2) to the muffin-tins of the two atoms 
participating in the covalent bonding. The FWF2s were 
nearly identical with the MLWFs, having the same cen- 
ters and negligibly different spreads, in particular. The 
FWFls are not centered in the middle of the C-C-bond, 
the FWF2s are, however, centered. Irrespective of the 
starting point (i.e. either FWF1 or FWF2) we arrive at 
the same MLWFs, which are bond centered. 

Figure [5] shows the contour plot of one of the three sp 2 - 
bonds for the first-guess FWF1 and for the MLWF. Fig- 
ure [H shows the 7r-orbital. Centers and spreads are given 
in Table fVl The initial spread of 17.08 a.u. 2 characteriz- 
ing the first-guess FWF1 is reduced by the minimization 
procedure to a final total spread of 16.23 a.u. 2 . 

The hopping matrix elements if m . m '(Ri — R2), 
Eq. (|2"9"]l . are listed in Table IVIl There is no coupling 



O (pz) 


3.730 


3.730 


0.629 


0.181 


O (px) 


3.730 


3.730 


0.686 


0.238 


O (py) 


3.730 


3.730 


0.686 


0.238 


O (pz) 


3.730 


0.000 


4.296 


0.012 


O (px) 


3.730 


0.000 


4.300 


0.016 


O (py) 


3.730 


0.000 


4.255 


-0.029 


O (pz) 


0.000 


3.730 


4.296 


0.012 


O (px) 


0.000 


3.730 


4.255 


-0.029 


O (py) 


0.000 


3.730 


4.300 


0.016 


Ba (pz) 


0.000 


0.000 


-0.047 


-0.047 


Ba (px) 


0.000 


0.000 


-0.011 


-0.011 


Ba (py) 


0.000 


0.000 


-0.011 


-0.011 


0(a) 


3.730 


3.730 


0.542 


0.095 


0(a) 


3.730 


0.000 


4.305 


0.021 


0(8) 


0.000 


3.730 


4.305 


0.021 


Ba (s) 


0.000 


0.000 


0.000 


0.000 


Ti (pz) 


3.730 


3.730 


3.863 


-0.038 


Ti (px) 


3.730 


3.730 


3.905 


0.003 


Ti (py) 


3.730 


3.730 


3.905 


0.003 
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FIG. 5: Contour plot of the FWFl (left) and MLWF (right) 
of an sp 2 -bond of graphene. 




FIG. 6: Isosurface plot of the 7r-orbital of graphene. 
Isosurface=±0. 1 

between the tt WFs and the sp 2 WFs. 



D. Platinum 

We close the results section with the discussion of the 
MLWFs for a Platinum chain. Our calculations were per- 
formed with the one-dimensional version^ of the FLEUR 
program and with spin-orbit couplin g 36 ' 37 i 38 ' 39 . The ex- 
tensions necessary to treat the spin-orbit case have been 
described in section III Gl The muffin-tin radii and the 
plane wave cut-off were chosen to be 2.22 a.u. and 3.7 
a.u. -1 , respectively. The RPBE^ exchange-correlation 
potential was used. The relaxed Pt-Pt distance is given 
by 4.48 a.u.. We calculated 12 MLWFs corresponding to 
the s- and d-states of Platinum using 8 fc-points. The 
localized trial orbitals were chosen to be eigenstates of 



TABLE V: Centers and spreads of the first-guess (first row) 
and maximally localized (second row) WFs (atomic units). 





X 


V 


z 


<x 2 > 


FWFl (sp 2 ) 
FWFl (sp 2 ) 
FWFl (sp 2 ) 
FWFl (tt) 


2.038 
2.038 
4.064 
2.714 


1.169 
-1.169 
0.000 
0.000 


0.000 
0.000 
0.000 
0.000 


2.184 
2.184 
2.184 
10.526 


MLWF (sp 2 ) 
MLWF (sp 2 ) 
MLWF (sp 2 ) 
MLWF (tt) 


2.035 
2.035 
4.070 
2.714 


1.175 
-1.175 
0.000 
0.000 


0.000 
0.000 
0.000 
0.000 


2.052 
2.052 
2.052 
10.075 



TABLE VI: Hopping matrix elements of graphene. Energies 
are in meV. 00, 10, 11 and 20 denote the translations of the 
obitals in units of the primitive translations. 





00 


10 


11 


20 


sp 2 (l),sp 2 (l) 
sp 2 (l),sp 2 (2) 
sp 2 (l),sp 2 (3) 
sp 2 (2),sp 2 (l) 
sp 2 (2),sp 2 (2) 
sp 2 (2),sp 2 (3) 
sp 2 (3),sp 2 (l) 
sp 2 (3),sp 2 (2) 
sp 2 (3),sp 2 (3) 


-15038 


560.7 


6.6 


51.3 


-2139 


78.0 


-21.5 


7.4 


-2139 


-144.1 


2.5 


-19.9 


-2139 


-529.8 


-21.5 


-21.5 


-15038 


-109.7 


6.6 


-6.7 


-2139 


78.0 


2.5 


7.4 


-2139 


-2139.1 


78.0 


-144.1 


-2139 


-529.8 


78.0 


-21.5 


-15038 


560.7 


-16.4 


51.3 


TT, TV 


-8329 


-728.0 


162.9 


51.6 



the z-projection of the spin operator. Both the direc- 
tion of the chain and the spin-quantization axis are given 
by the z-direction. We chose the angular parts of the 
trial-orbitals for the <i-bands to be d 3x 2_ r 2 1 d 3y 2_,,2, (i.e. 
rf 3z 2_ r 2 rotated to be coaxial with the x- and y-directions, 
respectively), d xy , d xz and d yz . The localized trial orbital 
corresponding to the sp-like WF was constructed as a lin- 
ear combination of two localized s-orbitals on neighboring 
atoms. The MLWFs are spinor-valued and complex. 6 
out of the 12 MLWFs are characterized by a dominance 
of the spin-up component while the spin-down compo- 
nent dominates the other 6 MLWFs. The two groups of 
spin-up and spin-down dominated WFs are symmetric 
by interchange of spins. Hence we will consider only the 
6 spin-up dominated WFs in the following, unless explic- 
itly stated. The angular dependencies of the real parts 
of the dominating spin-up components are approximately 
given by d xz and d yz , d% x 2_ r 2 and d 3y 2_ r 2, d xy , and sp. 
The MLWFs d xz , d yz and d 3x 2_ r 2, d 3y 2_ r 2 are symmetry 
equivalent, respectively. The sp-like WF is positioned 
bond-centred between two neighboring Pt-atoms. The 
angular functions that approximately describe the imagi- 
nary part of the spin-up component as well as the real and 
imaginary parts of the spin-down components agree very 
well qualitatively with the results of our simple model of 
section III Gl given in Table Q] We found qualitative de- 
viations only for the d 3y 2_ r 2 -orbital (and the symmetry- 
equivalent d 3a .2_ r 2-orbital) shown in Figure[7] While Ta- 
ble U predicts the real part of the spin-down component 
belonging to the d 3y 2_ r 2-orbital to vanish, it turns out to 
be non-vanishing and d^-like. This may be attributed 
to the fact that the actual G? 3l ,2_ r 2-like orbital is not rota- 
tionally invariant around the y-axis, but rather squeezed 
in x-direction. The d^-like WF is shown in Figure[8] As 
there is no spin-orbit coupling for s-states the spin-down 
component of the sp-like WF, which is shown in Figure 
[9l is p-like. 

Table IVLT1 lists the spreads. The maximal localization 
procedure reduces the initial total spread of 195.72 a.u. 2 
to a final total spread of 37.56 a.u. 2 . 
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FIG. 7: d 3a 2_ r 2-like orbital of a one-dimensional Pt-chain. 
Top row: Left: Real part of spin-up component (d 3y 2_ r 2, 
Isosurface=±0.1), Right: imaginary part of spin-up compo- 
nent (d X y, Isosurface=±0.001). Bottom row: Left: Real 
part of spin-down component (d xz , lsosurface=±0. 00073), 
Right: imaginary part of spin-down component (d yz , 
Isosurface=±0.0025). 

FIG. 8: d xy -\ike orbital of a one-dimensional Pt-chain. 
From left to right: Real part of spin-up component (d xy , 
Isosurface=±0.2), imaginary part of spin-up component 
(d x 2_ y 2, Isosurface=±0.005), real part of spin-down compo- 
nent (d yz , Isosurface=±0.001). 



In TableEEHwe list the spin-resolved nearest neighbor 
hopping matrix elements for the spin-up dominated ML- 
WFs between identical orbitals calculated according to 
Eq. ([55)1 . As the (|,|) components scale quadratically 
with the admixture of spin-down to the spin-up domi- 
nated WFs, they are small. Likewise the (T,l) compo- 
nents are found to be small: The angular distributions of 
the spin-down components of the WFs differ from those 
of the spin-up components, the admixture of spin-down is 



TABLE VII: Platinum chain: Spreads of the MLWFs (atomic 
units) . 





dxz 


d 3x 2_ r 2 


dxy 


sp 


<x 2 > 


3.336 


2.416 


2.326 


4.952 



TABLE VIII: Platinum chain: Spin-resolved nearest neighbor 
hopping matrix elements for the spin-up dominated MLWFs 
between identical orbitals (meV). 





dxzi dxz 


d 3x 2_ r 2, d 3x 2_ r 2 


dxy-, dxy 


s, s 


T,T 


1170.9 


-548.8 


-269.7 


-2481.7 


U 


-0.1 


0.4 


-0.1 


29.3 


•I) -1 


1.0 


-0.6 


-0.7 


-21.3 



small, and the spin-orbit coupling, which couples the two 
spin-channels, is important only close to the nuclear cores 
and hence the coupling between functions well-localized 
on different atoms is small. For the on-site hopping ma- 
trix elements, however, the (T,J.)- or (|, |)-components 
can dominate, because the two WFs are centered on the 
same atoms in this case, and their overlap close to the 
nuclear cores can be large. In Table IIXI we list a se- 
lection of spin-resolved on-site hopping matrix elements 
that are dominated by hopping from spin-up into spin- 
down, which is mediated by spin-orbit coupling. d xz is a 
spin-up dominated d^-like WF. According to Table|T]the 
spin-orbit interaction provides a coupling to d x 2_ y 2\ J,), 

which overlaps with dl x2 2 - Analogously, there is a 

transition from d^ y2 _ r2 to d yz \ |), which overlaps with 

d yz . The other two examples in Table IIXI are easily in- 
terpreted analogously on the basis of Table [J The (J,, f)- 
contributions in Table IIXI are negligibly small, because 
the | J,)- and | |)-components of the spin-up and spin- 
down dominated WFs are small, respectively. Table IXl is 
analogous to Table IVTlTl but now for the nearest neighbor 
hoppings. The comparison of the two Tables shows that 
the (|, |)-contributions decay fastest, which is consistent 
with the facts that the spin-orbit coupling is strongest 
close to the nucleii, and that the WFs are well localized. 




FIG. 9: sp-like orbital of a one-dimensional Pt-chain. Left: 
real part of spin-up component (sp, Isosurface=±0.04), Right: 
real part of spin-down component (p x , Isosurface=±0.004). 



TABLE IX: Platinum chain: Spin-resolved on-site hopping 
matrix elements between spin-up and spin-down dominated 
MLWFs (meV). 





d} d l 

U X Z, u Sx 2_ r 2 


(7 T d l 

LL - iv 2_ r 2 , u yz 


d} d l 

x z ) ^xy 


d T d l 

^xy riz 


T,T 


-142 


134 


10 


-6 


U 


460 


460 


268 


268 


l,T 














-1 1 -1 


134 


-142 


-6 


10 
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TABLE X: Platinum chain: Spin-resolved nearest neigh- 
bor hopping matrix elements between spin-up and spin-down 
dominated MLWFs (meV) . 





U X2) U, % x 2_ r 2 


u 3l) 2_ r 2 , u yz 


^ajz) xy 




u 


33 


0.8 


5.6 


-9.0 


u 


9.8 


9.8 


7.5 


7.5 


u 














•i ) i- 


0.8 


33 


-9.0 


5.6 



IV. CONCLUSIONS 

We have described the implementation of Wannier 
functions within the FLAPW program FLEUR for bulk, 
film and wire geometry. Two kinds of WFs with op- 
timized localization properties - the first-guess and the 
maximally localized Wannier functions - have been de- 
scribed and calculated for four concrete systems, SrV0 3 , 
BaTiOa, graphene and platinum. Our results are in very 
good agreement to previous ones, where available, includ- 
ing the ferroelectric polarization of BaTiC^. We found 
the first-guess WFs and the MLWFs to be similar for the 
first three systems, and rather different for Pt. While in 
cases where the first-guess WFs and the MLWFs do not 
differ substantially there is the option to use the first- 
guess WFs in practice for certain applications, which is 
computationally less demanding, the extended scheme 
needed for the construction of the MLWFs still proves 
valuable if quantities such as the electric polarization are 
supposed to be extracted. 
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APPENDIX A: 

THE Mi?» 



VACUUM CONTRIBUTIONS TO 
MATRIX IN CASE OF FILM 
CALCULATIONS 



In case of the film implementation of the FLAPW 
method, an additional semi-infinite vacuum region is 
present, which results in an additional contribution to 
the wave function overlaps Mmn\vAC- In this appendix 
we give explicit expressions for the vacuum contributions 
to the Mmn 3 ^ matrix elements. 

In the film geometry, the interstitial region stretches 
in z-direction from —D/2 to D/2, which is chosen to be 
the direction orthogonal to the film. Thus, one of the 
two vacua extends from — oo to —D/2 while the second 



vacuum extends from D/2 to +oo. The two vacua are 
treated analogously and we will restrict the discussion 
to the vacuum between D/2 and +oo. According to the 
topology of the vacuum region, the Bloch wave functions 
in the vacuum are represented in the following way: 



^k|,m(x)|vAC 



£*G„(k||.*)e <{G » +kl 



(Al) 



with 

(k„ , z) = (k,| (z) + (k||)u^ (z), (A2) 

where G = (Gm,G z ) and x = (xm,z) have been used, 
with G|| and xy the in-plane components. The fc-point k|| 
belongs to the two-dimensional BZ. ife (z) and tig (z) 
are the solution of the one-dimensional Schrodinger equa- 
tion in the vacuum and its energy derivative, respectively. 
Substituting Eq. |XT] into Eq. [T2] yields: 



M, 



(k,i,b) 



= J2 f „(k„,*))**2. ,([k | |+b],z)rf 3 x 

J—L, J VAC 11 

II ' I 

(A3) 

with Q = G|, - G|| - G(k|| + b). While vectors ky 
and [ky + b] always lie in the two-dimensional Brillouin 
zone, the b and G (ky + b) vectors have a z-component 
in general, which leads to the following expression for the 

Mmn' h ^ matrix elements: 



M (k„,b) 

lvl mn 



g,,,g: 



POO 

x / e- iG '< k ii +b >'(*S (k||, ; ([k t |+b],z)dz, 

Jd/2 11 

(A4) 

with 5|| being the in-plane unit-cell area, and the last 
integral is a linear combination of one-dimensional inte- 
grals of the form 



D/2 



D/2 



- iGa (k ll+ b)* n kj, ( z )u [ £} +h] (z)dz 



(A5) 



which are easily computed numerically for every pair of 
(Gii,G(|). 



APPENDIX B: VACUUM CONTRIBUTIONS TO 
THE M^ b) MATRIX IN CASE OF ONE 
DIMENSIONAL CALCULATIONS 

In the case of the one-dimensional setup the vacuum 
region surrounds a cylinder with the symmetry axis along 
the z-direction and radius i? vac . The wave function in the 
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vacuum is represented in the following form (in the ID 
case the Bloch vector is k = (0, 0, k z )): 

V*. m M = E(^;«^(A,,r) + B^ fc ;«^(fc»,r))x 

xe i PV e i(G i +fc i )z 

. ( B1 ) 

where x = (z,r,<p) in cylindrical coordinates, G z is the 
^-component of the reciprocal vector G, and p is an inte- 
ger number labeling a cylindrical angular harmonic. The 
exponentially decaying functions u and u are the solu- 
tions of the radial equation for the vacuum and its energy 
derivative, respectively. Taking into account the expan- 
sion of a plane wave in cylindrical coordinates 



e iGx = e iG ' z 



iPe^-^ljpiGrr), (B2) 



with S = 2ttT, and T standing for the lattice constant of 
the system under consideration along the z-axis. 



APPENDIX C: LOCAL ORBITAL 
CONTRIBUTIONS TO THE M^ b) MATRIX 



In order to increase the variational freedom of the 
FLAPW-basis or to describe semicore levels adequately, 
it may be supplemented by local orbitals4i In this case 
the expressions for the BFs in the spheres are modified: 



with tpQ and G r being cylindrical angular and radial co- 
ordinates, respectively, of the vector G = (G z ,G r ,fa) 
in reciprocal space, and J p standing for the cylindrical 
Bessel function of order p, the ID-vacuum contribution 
to the Mmn matrix reads: 



M (fe *' b) | 



-,„„ IVAC = / e 1 ' x (V^(x))>[fc*+b],n( x ) 

/VAC 



= E E 

G*,G> Z p,p> JVAC 

x e -tG||(k,+b)xn e i{v'- P )v *™ p ?£;(fc z , [k z +b],r) d 3 x, 

" (B3) 
where in analogy to the case of the film geometry, vectors 
b and G(k z + b) may have a non-zero component in 
the plane normal to the z-axis, and the function \P is 
constructed from the products of the u- and u-functions 
with corresponding A- and B-coefficients at fc-points k z 
and [k z + b}. Introducing the vector Q = G z — G z — 

G z (k z +h) the expression for the M^' b ' can be reduced 
to 



VAC 



E E- s '-^<- ' ' 



P' p - i (p-p')VG(k z 



x / r,J p ,_ p (G r (k z + b)r)*™£;gr (k z1 [k z + b], r) dr, 

J i?vac 

(B4) 



^ km (x)| MT , = ^(^, m (k)<(r) + J B^ m (k)<(r))F i (f) 
+ E^L, m (k)< Mr i0 (f), 

Lo 

(CI) 

where Lo = [lo, mo) stands for the corresponding values 
of the angular quantum numbers (I, m) assigned to each 
local orbital. Due to the local orbitals, additional terms 
arise in the expression Eq.QUfor theMi k ^ b) | 

mt' matrix: 

» / /-( k > b )|Lo _a -ib-T^v 

x ( E (^, m ( k ))*CLo', m ([k + b])^ (b, L, Lo')+ 

L,Lo' 

+ E (^, m (k))*CL,, m ([k + b])^(b,L,L ') + 

L,Lo' 

+ E (CL, m (k))*^, m ([k + b])^(b,L ,i') + 

Lo,L' 

+ E (CL,„(k))*^, m ([k + b])^ 2 (b,L ,L')+ 

Lo,L' 

+ E (CL, m (k))*^, m ([k + b])^(b,L ,L ')), 



Lo,Lo' 



(C2) 

where the corresponding radial function for the local or- 
bital is taken in the ^-integrals, whenever a radial func- 
tion u has an index lo. 
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